ESS 575 Models for Ecological Data

Dynamic Models: Forecasting Effects of Harvest on Lynx

April 24, 2019

---

Problem

The Eurasian lynx (Lynx lynx) is a medium-sized predator with broad distribution in the boreal forests of Europe and Siberia. The lynx is classified as a threatened species throughout much of its range. There is controversy about legal harvest of lynx in Sweden. Proponents of harvest argue that allowing hunting of lynx reduces illegal kill (poaching). Moreover, Sweden is committed to regulate lynx numbers to prevent excessive predation on reindeer because reindeer are critical to the livelihoods of indigenous pastoralists, the Sami. Many environmentalists oppose harvest, however, arguing that lynx are too rare to remove their fully protected status. A similar controversy surrounds management of wolves in the Western United States.

Your task is to help resolve this controversy by developing a forecasting model for the abundance of lynx in a management unit in Sweden. You have data on the number of lynx family groups censused in the unit as well as annual records of the number of lynx harvested from the unit. You will model the population using the deterministic model: \[N_t=\lambda(N_{t-1}-H_{t-1}).\] where \(N_{t}\) is the true, unobserved abundance of lynx and \(H_{t-1}\) is the number of lynx harvested during \(t-1\) to \(t\). The parentheses in this expression reflect the fact that harvest occurs immediately after census, such that the next years population increment comes from the post-harvest population size. (For the population ecologists: What would be the model if harvest occurred immediately before census? Three months after census? Continuously throughout the year?)

*Immediately before census:

\[N_t=\lambda N_{t-1}-H_{t-1}.\] *Three months after census:

\[N_t=(\lambda^{\frac{1}{4}}N_{t-1}-H_{t-1})\lambda^{\frac{3}{4}}\] \[N_t=(\lambda^{\frac{1}{4}}N_{t-1}-H_{t-1})\lambda^{\frac{3}{4}}\] *Throughout the year: \[N_t=(\lambda^{\frac{1}{2}}N_{t-1}-H_{t-1})\lambda^{\frac{1}{2}}\] \[N_t=\lambda N_{t-1}-\lambda^{\frac{1}{2}}H_{t-1}\]

We can reasonably assume the harvest is observed without error. It may be a bit of a stretch to assume an that observations of \(N_t\) are made perfectly, but my Swedish colleagues have convinced me that their census method, snow tracking, does a good job of estimating the number of family groups (if not the number of lynx) in a management region.

The challenge in this problem is that the observations of lynx abundance (family groups) is not the same as the observation of harvest (number of lynx). Fortunately, you have prior information on the proportional relationship between number of family groups and number of lynx in the population, i.e, \[\phi=\frac{f}{N}\] where \(f\) is the number of family groups and \(N\) is the population size, mean \(\phi\)=.163 with standard deviation of the mean = .012.

1. Develop a hierarchical Bayesian model (also called a state space model) of the lynx population in the management unit. Diagram the Bayesian network of knowns and unknowns and write out the posterior and factored joint distribution.
   

2. Write JAGS code to approximate the marginal posterior distribution of the unobserved, true state over time (\(\mathbf{N}\)), the parameters in the model \(\lambda\) and \(\phi\) as well as the process variance and observation variance. Summarize the marginal posterior distributions of the parameters and unobserved states. Data are in the accompanying file, Lynx data new.csv.

3. Check MCMC chains for model parameters, process variance, and latent states for convergence. This will probably require using the excl option in MCMCsummary.

4. Conduct posterior predictive checks by simulating a new dataset for family groups (\(f_t\)) at every MCMC iteration. Calculate a Bayesian p value using the sums of squared discrepancy between the observed and the predicted number of family groups based on observed and simulated data, \[T^{observed} = \sum_{t=1}^{n}(f_{t}^{observed}-N_{t}\phi)^{2} \\T^{model} = \sum_{t=1}^{n}(f_{t}^{simulated}-N_{t}\phi)^{2}.\] The Bayesian p value is the proportion of MCMC iterations for which \(T^{model}>T^{obs}\).

5. Assure that your process model adequately accounts for temporal autocorrelation in the residuals, allowing the assumption that they are independent and identically distributed. To do this, include a derived quantity \[e_t=y_t-N_t\] in your JAGS code and JAGS object. Use the following code or something like it to examine how autocorrelation in the residuals changes with time lag. Write a paragraph describing how to interpret the plot produced by this function.

acf(unlist(MCMCpstr(z,param="e",func=mean)),main="", lwd = 3, ci=0)

The autocorrelation of time series is the Pearson correlation between values of the process at two different times, as a function of the lag between the two times. It can take on values between -1 and 1. Values close to 1 or -1 indicate a high degree of correlation. The residuals not auto correlated if their values drop close to 0 at relatively short lags and alternate between positive and negative values at subsequent lags. This plot reveals no autocorrelation in the residuals of our model.

6. Plot the median of the marginal posterior distribution of the number of lynx family groups over time (1998-2016) including a highest posterior density interval. Include your forecast for 2017 (the predictive process distribution) in this plot. Your plot should resemble Figure 1, below.

7. Decisions on the number of lynx to be harvested must be made before the population is censused, even though harvest occurs immediately after census. This is because it takes months to properly issue licenses to the designated number of hunters. Make a forecast of the number of family groups in 2018 assuming five alternative levels for 2017 harvest (0, 10, 25, 50, and 75 animals). Environmentalists and hunters have agreed on a acceptable range for lynx abundance in the unit, 26 - 32 family groups. Compute the probability that the post-harvest number of family groups will be below, within, and above this range during 2018. Tabulate these values. Hint: Set up a “model experiment” in your JAGS code where you forecast the number of lynx family groups during 2018 under the specified levels of harvest. Extract the MCMC chains for the forecasted family groups( e.g., fg.hat) using MCMCchains(coda_object, params = fg.hat) Use the ecdf function on the R side to compute the probabilities that the forecasted number groups will be below, within, or above the acceptable range.

A note about the data. Each row in the data file gives the observed number of family groups for that year in column 2 and that year’s harvest in column 3. The harvest in each row influences the population size in the next row. So, for example, the 2016 harvest influences the 2017 population size.

Use a lognormal distribution to model the true lynx population size over time. Use a Poisson distribution for the data model relating the true, unobserved state (the total population size) to the observed data (number of family groups).

An alternative approach, which is slightly more difficult to code, is to model the process as \(\text{negative binomia}(N_t|\lambda(N_{t-1}-H_{t-1}, \rho))\) and model the data as \(\text{binomial}(y_t|N_t,\phi)\). Explain why this second formulation might be better than the formulation you are using. (It turns out they give virtually identical results.)

There are two advantages to the negative binomial process model and binomial data model. A negative binomial process model treats the true state as an integer, which for small populations like this one, has some advantages because it includes demographic stochasticity. The binomial data model assures that the observed state is never larger than the true state, which makes sense if the only source of error in the census is failing to observe family that are present. On the other hand, if there is a possibility of double-counting, which is the case here, the then Poisson is a better choice for the data model.

Code

Preliminaries

library(rjags)

## Loading required package: coda

## Linked to JAGS 4.2.0

## Loaded modules: basemod,bugs

library(MCMCvis)
library(HDInterval)
y=read.csv("Lynx data new.csv")
# Levels of  Harvest to evaluate relative to goals for forecasting part.
h=c(0, 10, 25, 50, 75)
#Function to get beta shape parameters from moments
shape_from_stats <- function(mu = mu.global, sigma = sigma.global){
         a <-(mu^2-mu^3-mu*sigma^2)/sigma^2
         b <- (mu-2*mu^2+mu^3-sigma^2+mu*sigma^2)/sigma^2
        shape_ps <- c(a,b)
        return(shape_ps)
}

#get parameters for distribution of population multiplier, 1/p
shapes=shape_from_stats(.163,.012)
#check prior on p using simulated data from beta distribution
x = seq(0,1,.001)
p=dbeta(x,shapes[1],shapes[2])
plot(x,p,typ="l",xlim=c(.1,.3))

Simulate data for initial values and model verification

I almost always simulate the true state by choosing some biologically reasonable values for model parameters and “eyeballing” the fit of the true state to the data. I then use these simulated values for initial conditions as you can see in the inits list below. This is particularly important. Failing to give reasonable initial conditions for dynamic models can cause no end of problems in model fitting. Remember, you must have initial conditions for all unobserved quantities in the posterior distribution. It is easy to forget this because some of these quantities don’t have priors.

##visually match simulated data with observations for initial conditions
#visually match simulated data with observations for initial conditions
endyr = nrow(y)
n=numeric(endyr+1)
mu=numeric(endyr+1) #use this for family groups
lambda=1.1
sigma.p=.00001
n[1] = y$census[1]

for(t in 2: (endyr+1)){
    n[t] <- lambda*(y$census[t-1] - .16 * y$harvest[t-1])  #use this for family groups
    }
plot(y$year, y$census,ylim=c(0,100),xlab="Year", ylab="Population size", main="Simulated data")
lines(y$year,n[1:length(y$year)])

Initial values and data

#Data for JAGS
data = list(
    y.endyr = endyr,
    y.a=shapes[1], 
    y.b=shapes[2],
    y.H=y$harvest,
    y=y$census,
    h=h
)

inits = list(
    list(
    lambda = 1.2,
    sigma.p = .01,
    N=n
    ),
    list(
    lambda = 1.01,
    sigma.p = .2,
    N=n*1.2),
    list(
    lambda = .95,
    sigma.p = .5,
    N=n*.5)
    )

{
sink("Lynx Harvest JAGS from cat.R")
cat("
model{
#priors==============
sigma.p ~ dunif(0,5)
tau.p <- 1/sigma.p^2
lambda ~ dunif(0,10)
p ~ dbeta(y.a,y.b)  #Get parameters a and b from mean and sd using moment matching to make this prior informative.  These are calcuated on R side and read in as data.

#Informative priors on initial conditions based on first year's observation of family groups
fg[1] ~ dpois(y[1])
N[1]~dlnorm(log(y[1]/p),tau.p)


#process model===============
for(t in 2:(y.endyr + 1)){  # the last year is a forecast with known harvest data
    mu[t] <- log(max(.0001,lambda*(N[t-1]-y.H[t-1])))
    N[t] ~ dlnorm(mu[t], tau.p)
    fg[t]<-N[t] * p
    }#end of process model
        
#data model===============
for(t in 2:y.endyr){   
         y[t] ~ dpois(p*N[t])  
    }  #end of data model

#Model checking
for(t in 1:y.endyr){
    #simulate new data for posterior predicitve check
        y.rep[t] ~ dpois(p*N[t])
        #accumlate test statistics for posterior predictive check
        sq[t] <- (y[t]-p*N[t])^2
        sq.rep[t] <-(y.rep[t] - p*N[t])^2
     #compute residuals for autocorrelation check
     e[t] <- y[t] - fg[t]
}
#calculate Bayesian P value
fit <- sum(sq[])
fit.new <- sum(sq.rep[])
pvalue <- step(fit.new - fit)

##forecast effects of different harvest regeimes on number of family grops during 2018
    for(i in 1:length(h)){
        #mu.hat is the forecast 1 year beyond y.endyr +1, i.e., 2019
        mu.hat[i] <- log(max(.001,lambda*(N[y.endyr+1]-h[i]))) 
        N.hat[i] ~ dlnorm(mu.hat[i], tau.p) #Nhat forecasts 2 years out
        fg.hat[i] <- N.hat[i] * p
    }
    
    
} #end of model

",fill=TRUE)
sink()
}

model = "Lynx Harvest JAGS from cat.R"

n.update=10000
n.iter=50000
n.adapt=5000
n.thin=1

jm = jags.model(model,data=data,inits=inits, n.adapt=n.adapt, n.chains=length(inits))

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 18
##    Unobserved stochastic nodes: 48
##    Total graph size: 1280
## 
## Initializing model

update(jm, n.iter=n.update)

z = coda.samples(jm,variable.names=c("lambda","sigma.p", "p", "e", "N", "N.hat", "fg", "fg.hat", "pvalue", "fit", "fit.new"), n.iter=n.iter, thin=n.thin)

#check convergence
MCMCsummary(z, excl = c("N.hat", "fg.hat", "fit.new", "fit", "e"), n.eff = TRUE)

##                mean          sd         2.5%         50%       97.5% Rhat
## N[1]    496.5156896 76.28075668 368.84907944 488.7399016 670.4303234    1
## N[2]    469.7096492 56.62811273 370.61732412 465.6803177 591.6167285    1
## N[3]    390.2884255 44.97743609 309.18146439 387.9006911 485.7912734    1
## N[4]    314.2946295 38.20409619 244.55029718 312.4519572 394.5381610    1
## N[5]    273.3476645 34.35699348 210.72687808 271.6233122 345.2568333    1
## N[6]    244.0335397 31.90695608 185.58033734 242.5963903 311.1066209    1
## N[7]    229.5783486 31.13484414 171.95601724 228.2248684 294.1655647    1
## N[8]    247.9859219 32.27583084 189.86465462 246.2948104 316.5002178    1
## N[9]    266.7357567 33.40998081 206.73022796 264.9062291 337.9481093    1
## N[10]   274.2936560 33.83509177 212.91235991 272.5158964 345.5465383    1
## N[11]   303.6993621 37.01857389 239.33360948 300.8676864 384.7242275    1
## N[12]   299.2845205 35.63054689 237.25337519 296.4002943 377.2784635    1
## N[13]   240.9234947 29.67477153 187.39244454 239.3626791 304.0042979    1
## N[14]   219.9985603 27.53173809 171.96059720 217.7507214 280.0927634    1
## N[15]   164.9039961 23.95170670 123.05034188 163.1754603 216.7133131    1
## N[16]   135.9671482 22.02427244  95.06800763 135.2165481 181.7185322    1
## N[17]   145.9454469 22.33655524 104.84895576 144.9705335 192.1275888    1
## N[18]   184.1861496 26.29410148 138.34998938 182.2486138 241.3034892    1
## N[19]   206.7630533 33.63033828 149.16478285 203.9490838 280.6228429    1
## N[20]   182.5466836 52.83240181 104.58150879 174.4802479 307.5605546    1
## fg[1]    81.9977600  9.06080326  65.00000000  82.0000000 100.0000000    1
## fg[2]    76.8476326  7.22099592  64.14476949  76.3308377  92.3027422    1
## fg[3]    63.8681965  5.72431941  53.23683152  63.6748917  75.7595328    1
## fg[4]    51.4287531  4.98084031  41.90668624  51.3460807  61.4975022    1
## fg[5]    44.7280515  4.55792109  36.07504406  44.6356716  53.9546399    1
## fg[6]    39.9304539  4.30616259  31.71147244  39.8452060  48.6076612    1
## fg[7]    37.5686553  4.29810753  29.22880178  37.5407601  46.0861438    1
## fg[8]    40.5749081  4.32097230  32.44794775  40.4491057  49.4513379    1
## fg[9]    43.6504901  4.46161889  35.28645930  43.5002559  52.8866212    1
## fg[10]   44.8926081  4.53624303  36.23920989  44.7921278  54.1958337    1
## fg[11]   49.7043191  4.92047890  40.93063598  49.3711146  60.3158673    1
## fg[12]   48.9898306  4.76216713  40.49760525  48.6644076  59.2571390    1
## fg[13]   39.4356611  4.02080213  31.80396849  39.3337776  47.7006955    1
## fg[14]   36.0066996  3.71411556  29.37728398  35.7640267  44.0719935    1
## fg[15]   26.9703974  3.26685284  21.09220364  26.7751973  33.9458066    1
## fg[16]   22.2444576  3.17208450  15.99416706  22.2572414  28.4621785    1
## fg[17]   23.8790437  3.18171994  17.71676088  23.8458351  30.2670714    1
## fg[18]   30.1408194  3.70395739  23.53445312  29.9089235  38.1060880    1
## fg[19]   33.8388961  4.94204786  25.26298286  33.4630952  44.5404092    1
## fg[20]   29.8579135  8.28847534  17.57330537  28.6170787  49.5382108    1
## lambda    1.0634226  0.04767295   0.97337322   1.0612756   1.1655843    1
## p         0.1645009  0.01194655   0.14189762   0.1641872   0.1886823    1
## pvalue    0.5777400  0.49392117   0.00000000   1.0000000   1.0000000    1
## sigma.p   0.1671483  0.06110365   0.07598267   0.1579512   0.3107395    1
##          n.eff
## N[1]     13653
## N[2]      9241
## N[3]      9476
## N[4]      9405
## N[5]      9515
## N[6]      9611
## N[7]      9879
## N[8]      9950
## N[9]     10669
## N[10]    10514
## N[11]    10216
## N[12]    10422
## N[13]    10332
## N[14]     9781
## N[15]     9783
## N[16]    10154
## N[17]    10171
## N[18]    11768
## N[19]    13915
## N[20]    26533
## fg[1]   150530
## fg[2]    32104
## fg[3]    62050
## fg[4]    47505
## fg[5]    42220
## fg[6]    36178
## fg[7]    26137
## fg[8]    47506
## fg[9]    55779
## fg[10]   56821
## fg[11]   50130
## fg[12]   47255
## fg[13]   48963
## fg[14]   42991
## fg[15]   34832
## fg[16]   18749
## fg[17]   24904
## fg[18]   33032
## fg[19]   29216
## fg[20]   43543
## lambda   42051
## p         5193
## pvalue   99297
## sigma.p   9458

#Do goodness of fit plot.
par(mfrow=c(1,1))
fit.new = MCMCchains(z,"fit.new")
fit =MCMCchains(z, "fit")
plot(fit.new,fit, xlab = "Discrepancy observed", ylab= "Discrepancy simulated", cex=.05, xlim=c(0,3000), ylim=c(0,3000))
abline(0,1)
p=MCMCpstr(z, params = "pvalue", func = mean)
text(500,2500, paste("Bayesian P value = ", as.character(signif(p$pvalue,2))))

acf(unlist(MCMCpstr(z,param="e",func=mean)),main="Autocorrelation", lwd = 3, ci=0)

#Plot true, unobserved state and forecast vs observations.
par(mfrow=c(1,1))
# +1 gives the one year forecast
years=seq(1998,y[nrow(y),1]+1)
fg = MCMCpstr(z, params = "fg", func=median)
fgHPDI = MCMCpstr(z, params = "fg", func = function(x) hdi(x, .95))
y2=c(y$census, NA)
plot(years,y2, ylim=c(0,100), ylab="Number of Lynx Family Groups", xlab="Years")
lines(years,fg$fg)
lines(years,fgHPDI$fg[,2], lty="dashed")
lines(years,fgHPDI$fg[,1],lty="dashed")
abline(h=26,lty="dotted", col="red")
abline(h=32, lty="dotted", col="red")

Figure 1. Median population size of lynx (solid line) during 1997-2016 and forecasts for 2017 with 95% credible intervals (dashed lines). Red dotted lines give acceptable range of number of family groups determined in public input process.

# Levels of  Harvest to evaluate relative to goals for forecasting part.
h=c(0, 10, 25, 50, 75)
#Acceptable limits on poplation size, determined by public input process.
lower = 26
upper = 32
p.in = numeric(length(h))
p.over =numeric(length(h))
p.under = numeric(length(h))

#get chains for harvest experiment
fg.hat = MCMCchains(z,"fg.hat")

for(j in 1:length(h)){
    p1 = ecdf(fg.hat[,j])(upper)
    p.under[j] = ecdf(fg.hat[,j])(lower)
    p.in[j] = p1 - p.under[j]
    p.over[j] = 1-p1
}

#trim to reasonable signficiant digits
p.under=signif(p.under,2)
p.in=signif(p.in,2)
p.over=signif(p.over,2)

alt.table = as.data.frame(cbind(h,p.under,p.in,p.over))
names(alt.table)=c("Harvest", "P(under)", "P(in)", "P(over)")
alt.table

##   Harvest P(under) P(in) P(over)
## 1       0     0.31 0.270   0.430
## 2      10     0.39 0.250   0.360
## 3      25     0.51 0.220   0.270
## 4      50     0.68 0.150   0.160
## 5      75     0.81 0.094   0.094

#save(alt.table, file="Harvest table.Rdata")